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Abstract 

A matrix-based approach to numerical integration of the DGLAP evolution equations is pre- 
sented. The method arises naturally on discretisation of the Bjorken x variable, a necessary proce- 
dure for numerical integration. Owing to peculiar properties of the matrices involved, the resulting 
equations take on a particularly simple form and may be solved in closed analytical form in the 
variable t = ]n(ao/a). Such an approach affords parametrisation via data x bins, rather than 
fixed functional forms. Thus, with the aid of the full correlation matrix, appraisal of the be- 
haviour in different x regions is rendered more transparent and free of pollution from unphysical 
cross-correlations inherent to functional parametrisations. Computationally, the entire programme 
results in greater speed and stability; the matrix representation developed is extremely compact. 
Moreover, since the parameter dependence is linear, fitting is very stable and may be performed 
analytically in a single pass over the data values. 
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1. INTRODUCTION 



Since the pioneering work of Duke and Owens |L| in the early eighties there has been an 
enormous investment in the phenomenological study of parton densities. As is well-known, 
the quark and gluon distributions inside hadrons acquire a scale (Q 2 ) dependence via higher- 
order PQCD corrections. Account of such dependence is necessary, both for correct analysis 
of experimental data and for various questions of theoretical importance. The evolution of 
the parton densities with Q 2 is governed by the DGLAP equations 0, |||, f|], the integro- 
differential nature of which hampers their numerical solution. In the literature there exist 
numerous techniques: such as, the so-called brute-force method ||, the use of Laguerre 
polynomials || [7J and solution in Mellin moment space with subsequent inversion. Short- 
comings common to almost all are the computer time required and decreasing accuracy for 
x -> 0. 

The precise motivation for developing this approach was an attempt to avoid the strong 
interplay between parameters (and x regions) to which standard approaches are prone. For 
example, typically, in any standard approach the overall normalisation is a function of all 
available parameters and thus, in particular, the low-x asymptotic power is prisoner to any 
fluctuations that may occur in the high-x region. As will emerge, the natural solution to 
this problem, i.e., to use function values themselves (binned in x and Q 2 ) as the parameters, 
leads to a gross simplification of the integro-differential equations involved. Indeed, suitable 
discretisation of the x variable immediately allows exact integration of the differential equa- 
tions in Q 2 and, moreover, most calculations may be performed prior to actual evolution 
or data fitting. A further enormous simplification of the equations is obtained by a suitable 
choice of the x bins — the matrices governing evolution are then "banded" lower triangular, 
with two crucial effects: they commute (and may be expressed in terms of sums of the unit 
matrix and nil-potent matrices); moreover, multiplication of such matrices is an order n 2 
operation (as opposed to the usual n 3 ) and they may be stored in very compact form, i.e., 
storage requirements are order n (as opposed to the usual n 2 ). 

Before proceeding it should be pointed out that an approach related to that proposed 
here has been presented by Santorelli and Scrimieri ||. Rather than comment immediately 
on the similarities and differences, it is probably more convenient to highlight such in the 
following, as and when appropriate. 

This paper is structured as follows: the following section contains the basic derivation 
of the approach to evolution while the subsequent covers its extension to singlet densities 
and higher-orders; in section f| numerical results are presented and section |5| describes the 
application to data fitting. Finally, some conclusions and prospects for the technique are 
presented. The appendices contain a brief demonstration of the special properties of the 
matrices involved and details of kernel integration. 

2. OUTLINE OF THE METHOD 

In the following the general form of the evolution equations is first presented and sim- 
plifying substitutions of variables are then applied. Finally, the matrix equation is derived 
and solved formally. 
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2.1. The Evolution Equations 



In the simplest case, non-singlet (NS) to leading logarithmic approximation (LLA), the 
evolution equations for parton densities take on the form: 

where f(x,Q 2 ) represents some NS parton density and P the corresponding splitting func- 
tion. The convolution form does allow exact InQ 2 integration if transformed to Mellin 
moment space. However, this requires knowledge of the function over the entire x region. 
Thus, if data analysis is the main aim, it is desirable to remain in x space. 

2.2. Simplifying Substitution of Variables 

Let us start by applying two simplifying variable transformations. It is convenient to 
substitute the scale variable Q 2 with the variable t: 

( = l»^f, (2.2) 

where Ql is some starting scale for evolution (typically a very few GeV 2 ). The second 
substitution is for the Bjorken variables x and y: 

M = ln-, v = ln-. (2.3) 

x y 



With these two substitutions, eq. (|2.lQ simplifies to 

df(u,t) 



dt 



dvP(u-v)f(v,t), (2.4) 



where we have taken into account that to LLA a s is proportional to 1/ In Q 2 and various 
coefficients have been absorbed into the definition of P. 



2.3. The Matrix Equation and a First Solution 

To examine the numerical approach, consider performing the right-hand side integral via 
a naive trapezoidal rule over subintervals of size h: 

%^ = £jW*(f), (2-5) 

k=l 

where the following rather obvious definitions have been made: 

u k = kh, f k (t) = f{u k ,t), P mk = hP(u m - u k ), (2.6) 

the typical vanishing of f(x) at x — 1 has been exploited and a factor one half in the last 
term of the series has been omitted. By noting that the sum in eq. (|2.5| ) runs only up to 
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m, one sees that the matrix P m k is lower triangular. Thus, finally, writing the equation in 
matrix form, one obtains 

f(t) = Pf(t), (2.7) 

where the dot indicates a derivative w.r.t. t. 

One would now, at least in principle, need only diagonalise P, via a matrix D say: left 
multiplication by D^ 1 would then result in 

D~ l f{t) =D- 1 PDD~ l i{t). (2.8) 

Defining f = D _1 f and the diagonalised matrix Pq = D~ 1 PD, the final simple form would 
therefore be 

f(t) = P D f{t). (2.9) 
The exact solution could thus be written down directly: 

/m(t)=e^/ m (0), (2.10) 

where the 7 m would be just the eigenvalues of the matrix P. Transforming back to the 
original basis, one would then have 

f m (t)=J2 elkt d mkfk(0). (2.11) 

k 

This would then be an exact solution in t of the differential equation, i.e., only the x variable 
having been discretised and treated numerically. 

It turns out that the eigenvalues are very close to one another and so diagonalisation is 
very nearly singular. A possible approach might be to avoid degeneracy by suitable choice 
of binning and also to work in extended precision on the computer; however, as we shall 
now show, this apparent obstacle may be turned in to a virtue owing to the nature of the 
matrices generated. 



2.4. A Better Solution to the Matrix Equation 

The choice of equally spaced bins in u leads, in fact, to total degeneracy. This impedes 
matrix diagonalisation and thus the solution must be obtained differently However, in com- 
pensation, a natural simplification occurs: the elements Pkm only depend on the differences 
k — m (see appendix [B| for a proof). Thus, first of all, such matrices may be stored nu- 
merically as vectors, there being only n independent elements. More specifically, we may 
define 

Pkm Vk—mi (2-12) 

where pk is just a vector of length n. Moreover, matrix multiplication, inversion etc. only 
require 0(n 2 ) floating-point operations as compared with the 0(n 3 ) for general nxn matrix 
manipulation. The banded triangular nature of the matrix provides a further remarkable 
simplification: the set of all such matrices is abelian. This means that the NS equation may 
effectively be solved algebraically. Thus, P may be treated as a c-number in eq. ( |2.7| ); the 
solution is then trivial: 

f(t) = e pt f(0). (2.13) 
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There is a further useful property of the matrices appearing here: if the diagonal elements 
of such a matrix are zero then the matrix is nilpotent; i.e., if Aij = a^j is such an nxn 
matrix with ao = then A n = 0. Separating P according to P = p$l + P, commutation of 
the two terms being trivial, we may then write 

f(t) =e Pot e pt f(0), (2.14) 

where, owing to the nilpotency of P, the factor exp[Pt] is now given by a finite sum of 
power terms. It is easy to see that its calculation indeed only requires 0(n 3 ) multiplications 
and is thus equivalent to a single normal matrix product. Moreover, any function admitting 
a power-series expansion (e.g., sin, cos, log, square-root, etc.) is evaluated similarly. The n 
terms of any such sum, of course, need only be calculated once and then stored (as an nxn 
matrix), providing for very rapid evaluation for any t. 

3. EXTENSION TO SINGLET AND HIGHER-ORDER 

In all modern analyses the extension to both singlet and higher logarithmic accuracy is an 
absolute necessity. Thus, we now turn first to the extension to the singlet case and secondly 
to higher-order QCD corrections. 

3.1. The Singlet Equations 

In the singlet case the equations become a 2x2 system: 

E(t) = PssS(t) + P S9 g(t), 

g(t) = P, s E(t) + P 99 g(t), (3.1) 

where the singlet-quark density is E = Y2/Qf( x ) an d the P a b are individually of banded 
triangular form. If we rewrite eq. (|3.1|) as the outer product of 2x2 matrices and the banded 
triangular kernels, P a b, 

F(t)=VF(t), (3.2) 
one immediately sees that the form of the original solution still holds: namely, 

F(t) = exp[Pt]F(0). (3.3) 

Of course, the matrix V, not being of the banded triangular type, does not permit direct 
evaluation of the exponential via a finite polynomial. However, the 2x2 matrix elements 
are commuting matrices which may still then be treated as c-numbers; the solution is thus 
trivially obtained by diagonalisation of the 2x2 system: 

F(t) = D exp[T D t] D- 1 F(0), (3.4) 

where Vd — P>~ X VD is diagonal in the 2x2 space. Note also that the transformation matrix, 
D, may be chosen such that D^ 1 = D (and is real). The upper- left and lower-right diagonal 
blocks may now be treated separately along the lines of the NS case. Note that, for the 
DGLAP kernels calculated in QCD, the matrix V is not singular; in the case that it were, 
the solution would in fact be trivial. 
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Alternatively, we may exploit the simple nature of 2x2 matrices. Let us use a basis of 
Pauli matrices: 



r 3 



1 






1 



(3.5) 



Writing V = PqI + P3T3 + P+r + + P_t_ = Pol + P? where the unit matrix and P now lie in 
the 2x2 space and the P a lie in x-bin space, we obtain 



and thus 



(P 3 2 + P + P_)l = P 2 1, 



cosh(Pt)l + sinh(Pt) P/P 



(3.6) 



(3.7) 



Note that although Po and P remain matrices in x-bin space they are also still banded lower 
triangular and thus the hyperbolic functions and square-root may again be evaluated as 
finite polynomials. We note in passing that, while both forms are stable, the evaluation of 
eq. ( pT4| ) requires slightly less computing time than does ( |3.7| ). 



3.2. Higher Logarithmic Accuracy 

Extending to higher logarithmic accuracy presents no obstacle in the NS case, indeed, it 
is trivial to all orders in a s . The evolution equation now becomes 

N 

q(t) = ^(«(g 2 )) n pWq(t), (3.8) 

n=0 

where the sum runs up to N, the order of the approximation. Since the matrices P^ are 
still all of the banded triangular type (i.e., commuting), an exponential solution of the form 
of eq. ( |2.13| ) is still valid for any N: 

q(t) = e p(0) * e pW t 1 -"*) e p(2) • ■ • q(0). (3.9) 

From the above form, it is also clear that the corrections actually factorise order-by-order 
and thus it is even possible to examine their effects separately and in parallel in a single 
calculation. 

Turning next to the singlet case, we encounter the first and only serious complication: in 

N 

F(t) = ^?W e- n 'F(t), (3.10) 

n=0 

the 2x2 matrices, P'"', do not commute and thus naive exponentiation cannot furnish the 
required solution. Let us examine the system in the next-to-leading-logarithmic approxima- 
tion (NLLA), which we rewrite as follows: 

F(*) = [A + Be-*] F(t). (3.11) 

If D is the matrix that diagonalises A at the 2x2 level, with the diagonal A' D = D~ l AD 
and B' = D~ 1 BD, then defining the vector G(t) via 

F(t) = D exp [A D t + B' D {\ - e' 1 )] G(t) (3.12) 
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where B' D is the diagonal part of B' and transforming to the variable s = e 1 we obtain 
dG(s) 



ds 



exp{2 [A' 3 In s + B' 3 (s - 1)] a 3 } [B' + a+ + G ( s )- ( 3 - 13 ) 



Defining upper and lower elements of G by G±, the coupled first-order equations trans- 
form into a pair of decoupled second-order equations: 

d 2 G±{s) xn fA 3 n ,\dG±{s 



ds^ 



±2 + B'^j + B' + B'_ G ± {s). (3.14) 



The solution may then be obtained in a relatively straight-forward manner via power-series 
substitution. However, while the banded triangular nature of the matrices certainly facil- 
itates evaluation, the series is infinite here and thus must be truncated. Fortunately, the 
compact form of the matrices permits cheap computer storage of a large number of terms. 

We note here that the natural expansion is a power series in s and not t. Indeed, while 
for small t the two are essentially equivalent, for large values (where s — > 0) the expansion 
in t contains only a part of the series expressed in terms of s. In the following section we 
shall comment further on this point with regard to the approach of ||. 



4. NUMERICAL RESULTS ON EVOLUTION 

We have thoroughly tested the approach to LLA for the both the NS and full singlet 
evolution equations: the performance is highly satisfactory both in terms of precision and 
computing time. For the purposes of testing we have used the following representative set 
of input distributions at Ql = 4GeV 2 : 

qNs(x,Q 2 ) = A NS x~ - 5 (1- x) 3 , 
Z(x,Qt) = Asx^il-x) 3 , (4.1) 
g(x,Ql) = A g x- l {l-xf*. 

The normalisations are fixed by various sum rules, which also provide means of cross-checking 
the accuracy of the method. 

As a first examination, let us consider NS evolution using just 20 points per decade (i.e., 
there are 20 bins between any x and x/10) and a simple two-point interpolation to evolve 
over a range of t — to 2; this corresponds to an astronomical final Q 2 = 2.3 x 10 12 GeV 2 
(taking Qq = 4 GeV 2 and Aqcd = 0.240 GeV for three active flavours) but we feel that 
this is important to test stability of the algorithm (we shall make further comments later). 
Under these extreme conditions (using double-precision arithmetic) the first moment of the 
valence quark distribution remains constant to one part in 10 5 while the second moments 
(the momentum sum rule) agree with the analytical evaluation to 1% (valence and singlet 
quark) and 2% (gluon). Note, however, that to obtain 10~ 5 precision in the first-moment 
integral it is necessary (for this number of bins) to cover the region x G (10~ 10 , 1) and to 
employ higher-order interpolation formulae for the moment integral itself. The precision 
degreases rapidly with increasing moment (owing to heavier sampling of the large- a; region) : 
it is down to the 10% level by the sixth (fourth) moment for the valence and singlet quark 
(gluon). 

A significant improvement is obtained by also using higher-order interpolation formulae to 
calculate the matrices. Note that, in order to maintain the matrix properties, it is necessary 
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to interpolate in each x bin using only function values at the bin boundaries and higher 
values of x. This means that fake zero values must be introduced for a very few x values 
larger than one (this essentially cancels the gain for the last bins) . Eight-point interpolation 
leads to a gain in precision of around three orders of magnitude for the second moments of 
all the distributions. 

Clearly, a substantial gain may also be made by increasing the density of the x points, 
this also permits a check of precision as a function of x; moments are only sensitive to a 
very limited range of x, around the peak of x n ~ l q(x). Doubling the density of points to 
40 per decade results in an improvement in precision of a factor 3 (50) for the two-point 
(8-point) interpolation. Doubling again to 80 points per decade (i.e., only 400 in the range 
x G (10 -5 , 1), the precision for 8-point interpolation is better than one part in 10 4 for all 
moments up to n = 12 and better than one in 10 6 for the momentum sum rule, the gain for 
the higher-order interpolation is still overwhelming. As a final step we double once more to 
160 points per decade, but find no further significant improvement. 

Thus, we choose to compare results with a benchmark of 80 points per decade using 
8-point interpolation. For evolution with 20 points per decade and 8-point interpolation, 
the precision from x = 0.25 down to x = 10~ 10 is better than 5 parts in 10 6 , for the valence 
density, still 2% for the highest value available for this binning (x = 0.8) but already 2.4 
in 10 4 for x = 0.5. The singlet-quark (gluon) density has a roughly factor 20 (30) poorer 
accuracy except for the higher x values where the singlet-quark precision is comparable 
to that of the valence. As a final test, we compare evolution with 10 points. Even here 
performance is excellent: the valence precision is better than 1 in 10 4 for x < 0.2 while that 
of the singlet quark and gluon is better than one per mille over most of the range. Note that 
if better precision in the high x region is necessary, then a separate (parallel) evaluation 
may be made with a greater concentration of points over a limited region (x > 0.2 say). 

With regard to the large value of t chosen for testing, we now comment on the method 
proposed in M: there the truncation of the power series in t at no more than 12 terms is 
justified on the basis of rapid convergence. Direct examination of the full finite expansion 
calculated here suggests this might be illusory: some coefficients do not become truly small 
until much later in the series and thus it may be that convergence is only achieved in || 
thanks to the not large values of t considered there. We have studied the problem only 
superficially, but find that truncation at small order of the series obtained here renders the 
procedure unstable and, moreover, the onset of instability is very sudden with increasing t. 
We should point out, however, the authors of |§ claim to see no evidence for such behaviour 
in their studies |§. 

It might be held that the typical size of t corresponding to physically accessible energy 
scales is not normally larger than about 0.3, as considered in ||. On the other hand, if 
one wishes to study the dynamical generation of parton distributions starting at some very 
small scale (see, for example, \Wl), then the range of t may be considerably larger. Finally, 
as commented above, a more suitable variable would be s; in this case the starting point 
being s = 1 requires expansion in the variable s — 1, which is still always a small variable 
(0<s-l<l). 

We have not yet implemented higher-order corrections; however, in the NS case the algo- 
rithm will suffer no change, indeed, examination of eqs. ( |3.12| , |3.13| ) reveals that higher-order 
corrections are implemented as further multiplicative (smaller) corrections. A possible strat- 
egy in the singlet case might be to exponentiate the leading corrections (possibly together 
with all diagonal pieces) and then apply the approach of M to the residual equations. Al- 



ternatively, as shown in Section |37|, after all possible simple factorisation the remaining 
differential equation is amenable to power-series solution. 



DATA FITTING 



It is now relatively straight-forward to perform data fitting. The / m (0) in eq. Q2.14 ) may 



be taken as the input parameters to fit, directly providing the parton density at some starting 
Qq. The exact nature of the solution in t means that to fit a data point at some given t it is 
possible to obtain the necessary f m at precisely that value and then interpolate in x to the 
required x value. Thus, since the matrix P is calculated once-and-for-all at the start of the 
fit (and only once even for many fits), the number of operations necessary to evaluate the 
fit function for a given set of values of x and Q 2 is drastically reduced with respect to all 
standard approaches. The number of floating-point operations may be reduced further by 
pre-calculating the evolution operator, e Pot e pt , for a set of t values and then interpolating 
both in x and t. Note that the compact matrix storage requires just a single vector of length 
n for each value of t. Note too that evolution may be performed using a larger number of x 
bins than are used for parametrisation. 

In a fit, since each evolved value, f(x,t), depends only linearly on the input starting 
values, / m (0) = f(x m ,0), x 2 is quadratic in the parameters, precisely the / m (0). Thus, 
minimisation and calculation of the full error matrix can, in principle, be performed with 
a single pass over the data points; this affords considerable saving in computer time and 
storage space. Indeed, data values need never be stored, the only memory requirements 
being for the evolution matrices (in vector form) and parameter vectors, together with the 
resulting covariance matrix. 

Moreover, fixing the overall normalisation (i.e., imposing sum rules) at this stage is en- 
tirely superfluous and so has no influence, while the fit will naturally return a full covariance 
matrix for the parameters f m (0), from which it is then possible to deduce whatever is of in- 
terest as to asymptotic behaviour, satisfaction of sum rules, etc., in an entirely independent 
manner. 

So far we have only tested a NS code by generating fake data (with nominal but realistic 
errors) at some fixed large Q 2 . The data generated were not smeared statistically so that any 
resulting non-zero x 2 would be due purely to inaccuracy in the evolution (based on a limited 
number of parameters). As an example we have used 200 data points generated over the 
interval x G [10 -5 , 1] and for the same extreme Q 2 as above, with 5 equally spaced parameter 
values per decade (note, however, that there is no restriction on the parameter spacing) and 
evolution using 20 points per decade. We find a total resulting x 2 — 0.6 (i.e., 0.004/DoF), 
the highest x parameter (corresponding to x — 0.63) is returned to within 8% of its true 
value while for all x < 0.1 the precision is better than 2 parts per mille. The entire fitting 
procedure (including data reading) takes around one second on a Pentium-driven laptop. 
As a final test, we generate an accurate set of 400 data points over four Q 2 slices (equally 
spaced in t) in the same x range and fit using the same number of parameters and evolution 
points as above. The results are very similar to the previous. 



9 



6. FINAL COMMENTS 



There remain certain technical issues to address, in particular, the matrix manipulations 
involved: at any point instabilities may, in principle, arise. However, this is unlikely; for h 
small the matrix solution will always be equivalent to other methods (and, in particular, it 
must be precisely equivalent to the so-called brute-force method) and will therefore have a 
well-defined and unique solution. Moreover, much of the procedure may be performed to 
essentially arbitrary precision beforehand, the resulting matrices then being used as fixed 
input to the evolution/fitting procedure. 

A further technical issue concerns the discretisation in x: the number of points will 
always be small, compared with what might normally be adopted for accurate numerical 
integration. However, the real data do actually lie in a finite number of finite-size bins and 
thus discretisation is in any case effectively forced even in methods where some functional 
form is fit to the binned data. 

Finally, as discussed, while the extension to NLLA in the NS case proves entirely trivial 
and to singlet in LLA is relatively straight-forward (it becomes a 2 x 2 system of equations), 
a full extension to the singlet case in NLLA proves less direct as the 2x2 system is non- 
commuting. However, the expansion that is to be truncated should be well-behaved insofar 
as the parameter is a s itself. 



7. NOTE ADDED TO COMPLETED MANUSCRIPT 

On completion of this manuscript the author became aware of unpublished work by 
Pascaud and Zomer [|11|, V^, which contains many of the basic ideas presented here. There 



are however certain comments we should like to make. Firstly, the discussion of the treatment 
of the high-x region is unnecessarily complicated; it is sufficient to simply choose a second x 
grid starting at x ~ 0.2 with n > 100 and not 1000 via the complicated formula suggested 
in [|12[| . Moreover, the interpolation adopted there is linear and no comment is made on 
the possibility of choosing higher-order polynomials, which, as noted here, do improve the 
precision considerably. Finally, and this is the most serious criticism: as discussed in the 



present work, a solution of the form of eq. (9) in [|l!| cannot be trivially extended to the 
2x2 set of singlet equations; the LLA and NLLA matrices do not commute, thus direct 
exponentiation is not possible. In other words, the expression for e Ak , despite having the 
simplified form given in the cited papers, does not correspond to the correct evolution 
operator for the singlet densities. Indeed, the authors present only limited cross-checking of 
their results using the analytically calculable moments of the parton distributions, relying 
purely on the assumed accuracy of the largest-ra calculation. 



APPENDIX A: THE DGLAP KERNELS IN NUMERICAL FORM 



The kernel matrices may be calculated by hand with no great difficulty. For simplicity, 
only the NS case will be treated here (the singlet case is more involved numerically but not 
conceptually). Thus, the equation to solve implies a convolution integral of some NS parton 
distribution, which is taken to be of the form[]l3| f(x) = xq(x), with the following kernel: 



P(x) oc x 



[l + x 1 



X) 



r 



(Al) 
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where the '+' regularisation is defined by 



dx ■ 



/(*) 



'1 



x 



(A2) 



A simplification is obtained by observing that, under the integral, the plus regularisation 
obeys the following identity: 



dy 



V 



1-3/ 



f{y) 



P 1 72 "I 

/ d l / 7 -?4-/(y) + (A3) 



The integral to evaluate may therefore be cast into the form 



dy 

y 



1+3/' 



dy 



dy 



1-3/ 
1+3/ 2 
1-3/ 
1 + 3/ 2 



/ 



1 



y 



A- 



/(- 



dy 



i + r 



2/ 

dy 



/(*) 



2/ 



/(*) 



(A4) 



We choose to compute the integrals using the x variable rather than the transformed variable 
u introduced in the main text for simplicity of calculation; note that since the rc-step ratio 
is always rather close to unity, there is negligible numerical difference between the two 
approaches. 

The integrals may be simplified via integration by parts to obtain 



2ln(l-y)+y + -y' 



f 



+ 



dy ( 21n(l - y) + y + -y 1 



dy 



f 



+ 



2ln{l -y)+y + -y 2 



in which careful examination reveals that the first and third terms completely cancel. Upon 
substitution of y — > x/y we finally obtain 



1 dy x 

y y 



In 1 - 



x 



x lx 

+ - + -- 

y %y 



dy 



/(!/)• 



(A6) 



The integral to be performed requires the division of the interval [x m i n , 1] into n sub- 
intervals of equal size in the variable u, introduced earlier. While one might directly apply 
a trapezoidal (or better) rule, since the integrand is a product of the unknown distribution 
function and a known splitting function (with singularities), it is clearly better to use the 
latter as a weight function. Thus, if we define the k-th interval as [xk,Xk-i], with x^ = X k 
(and x min = A"), then from eq. ( |2.4j) we have 



— - ' y y 



k=l 



ln( + ^ + 

y 2 y 1 



[fk ~ Jfc-l) 
(X k - Xk-l) 



(A7) 



11 



where we have identified x with x m and assumed for simplicity a two-point interpolation 
formula for f(y). The coefficients of fk for each x m are then just the elements P m k of the 
DGLAP kernel matrix required. 

It is immediately clear that precision will be poorer in the region x — > 1. However, two 
points should be borne in mind: firstly, the data are sparse in this region and thus high 
precision is superfluous, and secondly, the absolute value of q(x) there is so small as to be 
negligible. Moreover, as and when greater precision for x — > 1 might be necessary, at no 
more than the cost of doubling computer time, evolution could be performed in parallel for 
the entire region of interest and for just the region x G [0.2, 1], say. 

In concluding this appendix let us remark that the form in which the elements are cal- 
culated numerically is crucial for the avoidance of large round-off errors due to repeated 
division and multiplication by 1 — A, where A is defined just after eq. fl2.10|) . We note also in 
passing that as the interval size tends to zero, the coefficients diverge only logarithmically, 
as all pole terms explicitly cancel. We further note that the use of higher-order interpola- 
tion formula? leads to significantly greater precision: we have tested Lagrange polynomials 
up to order seven and found improvements in precision of rather more than two orders of 
magnitude. 



APPENDIX B: PROOF OF BANDED FORM OF THE DGLAP MATRICES 

Here we prove the statement that the matrix elements Pk m only depends on k — m. First 
of all we note that the plus regularisation involves the subtraction of an integral over the 
entire interval [0, 1], which therefore depends neither on m nor k, indeed, it contains no x 
dependence, except in the factor f(x) itself. Thus, for the proof we may ignore the subtleties 
of the plus regularisation, then the typical integral to evaluate takes the form 

ff-'^P^W (B!) 



k=l Jx k 



y v y 



The two-point interpolating (Lagrange) function for f(y) between the points x = Xk and 
x fc _i is just 

/(») « V ~* k fk-i + V ~_ Xk -' fk (B2) 
Making the substitution y = Xk-\£, and using x k = X h we obtain 



£/ f P( ^ ) /(A^O, (B3) 



k=l 

where the interpolation formula now becomes 

/(A*-ifl » f k _, + Lzl f k . (B4) 

Thus the final result for the coefficients of the fk, in the integral governing the evolution at 
x m , clearly only depend on the difference m — k, which in turn only appears in the numerator 
of the argument of P in eq. (|B3|) . Clearly, increasing the number of interpolation points 
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does not affect the argument, provided the set of points used for any given bin terminates at 
the lower x boundary of that bin. 
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